Decoding Chaos: Geospatial Data Exploratory Analysis (GDEA)
Author
Teo Suan Ern
Published
February 27, 2024
Modified
February 29, 2024
1. Overview
Note: This is a continuation of the take-home exercise 4. For Project Brief, Project Objective, Storyboard and Initial Exploratory Data Analysis, please view it under Initial Data Exploratory Analysis segment.
Use st_crs() to check the set coordinate system of myan_bound.
st_crs(myan_bound)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Use st_transform() to assign correct EPSG code to dataframe, from one coordinate system to another coordinate system mathematically.
myan_bound4144 <-st_transform(myan_bound, 4144)
st_crs(myan_bound4144)
Coordinate Reference System:
User input: EPSG:4144
wkt:
GEOGCRS["Kalianpur 1937",
DATUM["Kalianpur 1937",
ELLIPSOID["Everest 1830 (1937 Adjustment)",6377276.345,300.8017,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Geodesy."],
AREA["Bangladesh - onshore; India - mainland onshore; Myanmar - onshore and Moattama area offshore; Pakistan - onshore."],
BBOX[8.02,60.86,37.07,101.17]],
ID["EPSG",4144]]
Use st_join to create new variable (based on year) Use tapply to
To find out the number of fatalities and armed conflict incidents in each third-level administrative region.
We will extract the fatalities records by event_type field.
for (year in2010:2023) {# Filter data_sf4144 for the current year data_sf4144_year <- data_sf4144[data_sf4144$year == year, ]# Perform spatial join joined_data <-st_join(myan_bound4144, data_sf4144_year)############### Summarise total fatalities at each adm2 level in myan_bound4144# Replace NA with 0 for polygons with no fatalities total_fata_col <-paste0('total_fata_', year) myan_bound4144[[total_fata_col]] <-tapply(joined_data$fatalities, joined_data$shapeID, sum) myan_bound4144[[total_fata_col]][is.na(myan_bound4144[[total_fata_col]])] <-0# Summarise percentage of fatalities at each adm2 level in myan_bound4144 pct_fata_col <-paste0('pct_fata_', year) myan_bound4144[[pct_fata_col]] <-round(myan_bound4144[[total_fata_col]] /sum(myan_bound4144[[total_fata_col]]) *100, 2) myan_bound4144[[pct_fata_col]][is.na(myan_bound4144[[pct_fata_col]])] <-0############## # Summarise total incidents at each adm2 level in myan_bound4144# Replace NA with 0 for polygons with no incidents total_inci_col <-paste0('total_inci_', year) myan_bound4144[[total_inci_col]] <-tapply(joined_data$year, joined_data$shapeID, sum) myan_bound4144[[total_inci_col]][is.na(myan_bound4144[[total_inci_col]])] <-0# Summarise percentage of incidents at each adm2 level in myan_bound4144 pct_inci_col <-paste0('pct_fata_', year) myan_bound4144[[pct_inci_col]] <-round(myan_bound4144[[total_inci_col]] /sum(myan_bound4144[[total_inci_col]]) *100, 2) myan_bound4144[[pct_inci_col]][is.na(myan_bound4144[[pct_inci_col]])] <-0############## # Summarise total political violence rate at each adm2 level in myan_bound4144# Replace NA with 0 for polygons with no incidents total_politicalrate_col <-paste0('total_politicalrate_', year) myan_bound4144[[total_politicalrate_col]] <-tapply(joined_data$political_rate, joined_data$shapeID, sum) myan_bound4144[[total_politicalrate_col]][is.na(myan_bound4144[[total_politicalrate_col]])] <-0# Summarise percentage of political violence rate at each adm2 level in myan_bound4144 pct_politicalrate_col <-paste0('pct_politicalrate_', year) myan_bound4144[[pct_politicalrate_col]] <-round(myan_bound4144[[total_politicalrate_col]] /sum(myan_bound4144[[total_inci_col]]) *100, 2) myan_bound4144[[pct_politicalrate_col]][is.na(myan_bound4144[[pct_politicalrate_col]])] <-0############## # Summarise total civilian violence rate at each adm2 level in myan_bound4144# Replace NA with 0 for polygons with no incidents total_civilianrate_col <-paste0('total_civilianrate_', year) myan_bound4144[[total_civilianrate_col]] <-tapply(joined_data$civilian_rate, joined_data$shapeID, sum) myan_bound4144[[total_civilianrate_col]][is.na(myan_bound4144[[total_civilianrate_col]])] <-0# Summarise percentage of civilian violence rate at each adm2 level in myan_bound4144 pct_civilianrate_col <-paste0('pct_civilianrate_', year) myan_bound4144[[pct_civilianrate_col]] <-round(myan_bound4144[[total_civilianrate_col]] /sum(myan_bound4144[[total_inci_col]]) *100, 2) myan_bound4144[[pct_civilianrate_col]][is.na(myan_bound4144[[pct_civilianrate_col]])] <-0############## # Summarise non-state ex rate at each adm2 level in myan_bound4144# Replace NA with 0 for polygons with no incidents total_nonstateex_col <-paste0('total_nonstateex_', year) myan_bound4144[[total_nonstateex_col]] <-tapply(joined_data$non_state_exchange, joined_data$shapeID, sum) myan_bound4144[[total_nonstateex_col]][is.na(myan_bound4144[[total_civilianrate_col]])] <-0# Summarise percentage of non-state ex rate at each adm2 level in myan_bound4144 pct_nonstateex_col <-paste0('pct_nonstateex_', year) myan_bound4144[[pct_nonstateex_col]] <-round(myan_bound4144[[total_nonstateex_col]] /sum(myan_bound4144[[total_inci_col]]) *1000) myan_bound4144[[pct_nonstateex_col]][is.na(myan_bound4144[[pct_nonstateex_col]])] <-0############## # Summarise govt regain rate at each adm2 level in myan_bound4144# Replace NA with 0 for polygons with no incidents total_govtregain_col <-paste0('total_govtregain_', year) myan_bound4144[[total_govtregain_col]] <-tapply(joined_data$govt_regain_exchange, joined_data$shapeID, sum) myan_bound4144[[total_govtregain_col]][is.na(myan_bound4144[[total_govtregain_col]])] <-0# Summarise percentage of govt regain rate at each adm2 level in myan_bound4144 pct_govtregain_col <-paste0('pct_govtregain_', year) myan_bound4144[[pct_govtregain_col]] <-round(myan_bound4144[[total_govtregain_col]] /sum(myan_bound4144[[total_inci_col]]) *1000) myan_bound4144[[pct_govtregain_col]][is.na(myan_bound4144[[pct_govtregain_col]])] <-0############## # Summarise by event_type at each adm2 level in myan_bound4144 event_types <-unique(data_sf4144$event_type)for (event_type in event_types) { event_count_col <-paste0(event_type, '_', year) myan_bound4144[[event_count_col]] <-tapply(joined_data$event_type == event_type, joined_data$shapeID, sum) myan_bound4144[[event_count_col]][is.na(myan_bound4144[[event_count_col]])] <-0 }}